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We introduce a white graph expansion for the method of perturbative continuous unitary trans¬ 
formations when implemented as a linked cluster expansion. The essential idea behind an expansion 
in white graphs is to perform an optimized bookkeeping during the calculation by exploiting the 
model-independent effective Hamiltonian in second quantization and the associated inherent cluster 
additivity. This appoach is shown to be especially well suited for microscopic models with many 
coupling constants, since the total number of relevant graphs is drastically reduced. The white 
graph expansion is exemplified for a two-dimensional quantum spin model of coupled two-leg XXZ 
ladders. 
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I. INTRODUCTION 

High-order series expansions are an important tool 
in statistical physics. Typically, the linked clus¬ 
ter theorem is used to determine the correct ex¬ 
pression of physical quantities up to high orders 
in perturbation by performing calculations on finite 
linked clusters. Historically, such linked cluster ex¬ 
pansions date back to the 1950’s and 1960’s JTH4], 
where high-temperature series for extensive thermo¬ 
dynamic quantities have been determined. Similarly, 
extensive zero-temperature ground-state properties 
like the ground-state energy or susceptibilities can be 
calculated via linked cluster expansions as first pro¬ 
posed by Nickel in 1980 and implemented by Mar- 
land in 1981 [5]. Further progress has been made 
steadily over the years » 

For non-extensive quantities like excitation ener¬ 
gies, however, linked cluster expansions are more 
complicated and for a long time the calculation of an 
energy gap was only achieved by a deficient scheme 
which also needed unlinked clusters |7J. It took un¬ 
til 1996 when Gelfand set up a true linked cluster 
expansion for a one-particle dispersion |10| . Gelfand 
realized that quantum fluctuations of the vacuum 
have to be subtracted from one-particle hopping 
amplitudes to perform proper linked cluster expan¬ 
sions. Yet, unrealized at that time, this approach 
in terms of similarity transformations on graphs vi¬ 
olates the cluster additivity and is therefore inap¬ 
plicable when the ground state and the targeted 
excitation-subspace are characterized by identical 
quantum numbers. In 2000. it has been shown that 


the use of orthogonal transformations on graphs re¬ 
stores the cluster addiditivity and therefore allows 
for a consistent calculation of many-particle excita¬ 
tion energies mm- 

At the same time, an alternative route to linked- 
cluster expansions has been established by the 
method of perturbative continuous unitary transfor¬ 
mations (pCUTs) [14j allowing for the calculation of 
high-order series expansions for many-particle exci¬ 
tation energies as well as spectral densities |15j- I n 
contrast to the other approaches mentioned above, in 
pCUTs a quasi-particle conserving effective Hamil¬ 
tonian in second quantization is derived model- 
independently, with the constrained that the unper¬ 
turbed part of the Hamiltonian has an equidistant 
spectrum and is bounded from below. Interestingly, 
this method fulfills the cluster additivity by con¬ 
struction. In recent years, pCUTs were indeed used 
as linked-cluster expansions, i.e. a full graph decom¬ 
position has been implemented to calculate relevant 
matrix elements. Important examples are the deriva¬ 
tion of effective low-energy spin models for the Hub¬ 
bard model on the triangular and honeycomb lattice 
nsum, the calculation of the one-magnon gap for 
the fully-frustrated transverse field Ising model on 
the triangular and kagome lattice [18] , the treatment 
of topological phase transitions of perturbed non- 
Abelian string-net models on the honeycomb lattice 
mm, or the determination of the one-triplon dis¬ 
persion of coupled Heisenberg dimers on hypercubic 
lattices for arbitrary dimension m- 

Overall, linked cluster expansions constitute an ef¬ 
ficient tool with a vast variety of applications. At 
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this point, let us also mention the recently devel¬ 
oped non-perturbative variants of linked cluster ex¬ 
pansions combining a graph decomposition with ex¬ 
act diagonalization [221425] . with continuous unitary 
transformations m im m mi or density matrix 
renormalization group m- 

The numerical power of all (perturbative) linked- 
cluster approaches, especially in more than one di¬ 
mension, relies on a full graph decomposition. To 
this end the lattice is devided into small subclus¬ 
ters on which the actual calculations are carried out. 
Afterwards, results in the thermodynamic limit are 
obtained by embedding properly the finite-cluster re¬ 
sults into the infinite system. If all sites of two dif¬ 
ferent subclusters are coupled identically, these clus¬ 
ters are indistinguishable for the Hamiltonian and 
are called topologically equivalent. Therefore, the 
calculation can be restricted to topologically distinct 
clusters yielding a highly efficient approach. 

However, the number of topologically distinct 
graphs grows exponentially with the order of pertur¬ 
bation as well as with the total number of different 
coupling constants. Consequently, it becomes very 
hard to reach sufficiently high orders for problems 
with several expansion parameters. This is especially 
relevant for the comparison with experimental data 
where typically different coupling strengths are im¬ 
portant and have to be determined. 

In this work we introduce a new kind of graph 
expansion, a so-called white graph expansion , which 
overcomes the latter limitation to a great extent. 
The white graph expansion benefits directly from 
the underlying framework provided by the effective 
pCUT Hamiltonian in second quantization. 

The paper is organized as follows. We start by 
giving the set up of systems we focus on in this pa¬ 
per in Sect. |TT] In Sect. |III| we describe the pCUT 
method and its implementation as a linked cluster 
expansion. The white graph expansion is then intro¬ 
duced in Sect. [IV] and it is applied to coupled two-leg 
XXZ ladders in Sect.[V] Finally, we give conclusions 
in Sect. IVI1 


II. SET UP 


We consider a generic quantum lattice Hamilto¬ 
nian H at zero temperature. By decomposing the 
original lattice into a superlattice of supersites, one 


can always rewrite exactly Ti as 


N X 

n = n 0 + J2 x ^ U) • W 

3 = 1 

Here a supersite might be a spin, two linked spins 
like a dimer, or any other finite set of linked sites 
which can be easily diagonalized and which has an 
equidistant spectrum bounded from below. 

The unperturbed part of T-L is diagonal in super¬ 
sites i of the lattice and can be written as 

Ho = E 0 + J2fUi,a 

= E'o + Q , (2) 

where Eq denotes a constant and the sum over a 
runs over all excited local degrees of freedom. The 
unperturbed ground state |ref) is interpreted as the 
vacuum and is given as the product state |ref) = 
|0)... |0) with |0) being the lowest state of a super¬ 
site. Accordingly, fj a |ref) creates a local excitation 
of type a on supersite i. In the following we assume 
the local spectrum to be equidistant. Consequently, 
it is always possible to introduce the counting oper- 

ator Q = fii = Ei, a fi,Ji,cr 

Supersites interact via the perturbation 
v = Ej^ j) - The sum over j runs over all 
N\ perturbation parameters A j. The different oper¬ 
ators V ^ build the bonds of the lattice, so that one 
can assign a different n color" for each perturbation 
parameter A j. Here we restrict the discussion to 
Hamiltonians where the perturbation V couples 
two supersites. A generalization to perturbations 
which couple an arbitrary number of supersites is 
straightforward. 

As a consequence of Eq. ©, one can rewrite Eq.[l] 
as 


N 

H = H 0 + Y, Tn > ( 3 ) 

n=—N 

so that [Q,T n ] = nT n . Physically, the operator 
T n = JA A jTn^ corresponds to all operators where 
the change of energy quanta with respect to T~Lo is 
exactly n. The maximal (finite) change in energy 
quanta is called zb N. 
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III. PERTURBATIVE CONTINUOUS 
UNITARY TRANSFORMATIONS 

The pCUT method is an efficient tool for the cal¬ 
culation of high-order series expansions for Hamilto¬ 
nians of the form Eq. © [IMS]. The application 
consists essentially of two steps: i) the perturbative 
order-by-order solution of the so-called flow equa¬ 
tion yielding an effective Hamiltonian and effective 
observables and ii) a non-trivial extraction of the ef¬ 
fective low-energy physics. 

The first step depends only on the structure of the 
perturbation, namely, how many particles at most 
are created (annihilated) by the perturbation. In 
this sense, the obtained solution is model indepen¬ 
dent. However, this generality comes at the price of 
a complex model dependend extraction process be¬ 
cause the effective operators are not normal-ordered, 
i.e. the desired matrix elements have to be extracted 
in a separate (second) step. The real challenge 
lies in this second step which represents the model- 
dependend part and therefore typically is the bottle¬ 
neck of the calculation. The extraction process is 
based on the linked-cluster theorem which allows the 
determination of matrix elements in the thermody¬ 
namic limit by applying the effective operators to 
finite clusters. At the end of this process a normal- 
ordered low-energy Hamiltonian is obtained. 

Essentially, two different schemes for the last step 
are possible. First, in the standard scheme the cal¬ 
culation is carried out on a single cluster m M V 
[32] which is sufficiently large to contain all processes 
of a desired perturbative order. The implementa¬ 
tion of this approach is straightforward and details of 
the operator structure are irrelevant for the general 
scheme. However, in two or more dimensions, this 
is typically unfavourable because memory problems 
arise and a full cluster decomposition is the method 
of choice in. 

Let us also mention that a hybrid, a decomposi¬ 
tion into rectangular clusters, has been also succes- 
fully introduced recently for the pCUT method [33]. 
A rectangular cluster expansion seems to be espe¬ 
cially well suited for non-perturbative extensions of 
linked cluster expansions, since a typical length scale 
can be assigned to each rectangular graph, which is 
much less trivial when performing a full graph de¬ 
composition [23] [27]. 

In the following we start by briefly introducing the 
method of continuous unitary transformations and 


its perturbative variant, the pCUT method. After¬ 
wards, we formulate the pCUT method as a linked 
cluster expansion giving details to graph generation 
and graph embedding. 


A. Perturbative solution of the flow equation 

The objective of continuous unitary transforma¬ 
tions (CUTs) is to transform the Hamiltonian into an 
optimized basis representation [34] [35] • The Hamil¬ 
tonian is then considered as a continuous function 
of the flow parameter i with = 0) = H as 
the starting Hamiltonian and = oc) = H e ff as 
the effective Hamiltonian. Introducing the antiher¬ 
mit ian generator r](£) of the CUT, one gets the flow 
equation 

^p- = m,w)] ■ ( 4 ) 

Here we choose the quasi-particle generator m 

H3 136] 

= sgnfe - qj)H id {£) (5) 

in an eigenbasis | i) of the counting operator Q, 
i.e. Q\i) = qi\i), as introduced in Eq. ©• 

The commutator in the flow equation Eq. © leads 
typically to an infinite number of terms, so that 
a truncation must be performed. In the pCUT 
method, the truncation is done in a perturbative 
manner and, accordingly, the flow equation is solved 
via the perturbative ansatz 

H({\ j }-,£) = Ho+ ( 6 ) 

OO 

£ £ F(e-,m)T(m), 

£j k 3= k \m\=k 

implying, according to Eq. ©> 

OO 

^({A i };^)= £ A^.A^ £ 

£,• kj=k \m\=k 

F(£; m) sgn (M(m)) T(m). (7) 
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Here, the following notations have been introduced 


m = (mi, m 2 , 

,m k ) 

(8) 

nil G {0, ±1, d=2,.. 

,.,±N} 

(9) 

m = k 


(10) 

T(m) =T mi T m2 T m3 

■ ■ ■ Tm k 

(11) 

k 

M(m) = 


(12) 


Ansatz © leads to recursive differential equations 
for the real functions F(l\m ) with the initial con¬ 
ditions F(0,m) = 1 for |m| = 1 and F(0,m) = 0 
for |m| > 1 recovering the initial Hamiltonian. Due 
to the structure of the recursive equations, they can 
be solved in principal analytically. The calculation 
of the analytic functions F can be performed on a 
computer up to a certain order n max limited by the 
computation time and the memory usage only M- 
The solution leads to an effective Hamiltonian of 
the form 

oo 

n eS ({X j }) = Ho+ £ a £.. a ^ 

I Zj kj=k 

£ (13) 

\rn\=k rn i = 0 

with C(m ) G Q. The restriction JE m* = 0 reflects 
the particle-conserving property Q] = 0 of the 

final Hamiltonian. The summands can be viewed 
as virtual fluctuations of the new dressed particles 
defined by the effective Hamiltonian. 

In a similar manner one can calculate other physi¬ 
cal quantities with pCUT m ■ To this end the same 
unitary transformation has to be applied to any ob¬ 
servable O of interest 

= b? S ({Aj};^), (14) 

Using the perturbative ansatz 

oo 

0({Xj};£)= £ a £.. a ^ 

Ej kj=k 

k +1 

E E G(i\ m, i)0(m ; i) (15) 

i=l \m\=k 

with 

0(m;i) := T mi .. • E m?; _ x OT m . ... , (16) 



FIG. 1. All fluctuations involving the highlighted links 
are pooled together in this representation. The overall 
contribution of these fluctuations to the hopping element 
ti,j in the thermodynamic limit (left panel) can be identi¬ 
fied with a reduced contribution of a disconnected cluster 
(right panel) vanishing due to the linked cluster theorem. 


one obtains recursive differential equations for the 
functions G(^;m, i). The final result is given by 

oo 

o eS ({\j})= E a F--E? 

Ej kj=k 
/c+1 

E E G(m; i)0(m ; i), (17) 

i=l \m\=k 

with G(m, i) = G{£ = oo; m, i) G Q. In contrast to 
the effective Hamiltonian, effective observables are 
not quasi-particle conserving. 

The coefficients C{m ) and G(m; i) are indepen¬ 
dent of the model and can be straightforwardly ap¬ 
plied to all the problems matching the requirements 
of the pCUT method. This generality is only possi¬ 
ble at the expense of the second model-dependend 
extraction process described next for the effective 
Hamiltonian i7 e ff • The treatment of effective observ¬ 
ables can be done in complete analogy. 


B. pCUT as a linked cluster expansion 

Since the effective operators like H e ^ are not 
normal-ordered, matrix elements defining the effec¬ 
tive model are not directly accesible. In order to 
calculate the desired matrix elements, the effective 
Hamiltonian H e ff can be applied to a finite cluster 
yielding the correct contributions in the thermody¬ 
namic limit for a given order and matrix element. 
This crucial property results from the linked cluster 
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theorem which states that only linked processes have 
an overall contribution to cluster additive quantities. 


As outlined in Sect. [TTJ here we restrict the discus¬ 
sion to Hamiltonians where the perturbation couples 
always two supersites. The coupling between super¬ 
sites is represented by a link in the cluster and the 
operator can be decomposed as 


T n = Y Tn ’ 1 , (18) 

l 


where r n? / affects only the two supersites which are 
connected by the link l in the effective lattice. We 
say, that the operator r n? / acts on the link /. The 
sum runs over all links of the lattice. Inserting this 


decomposition into Eq. (13) yields 


(19) 


H eS ({\ j })=u 0 + y n n : 

Yj kj=k 

Y c ( m) Y k , 1 1 * 1~mk )lk 

\rn\=k Yi rrii=0 l\...lk 

In practice, however, the Hamiltonian is applied us¬ 
ing Eq. (13). Each summand of order k can be un¬ 


derstood as a virtual fluctuation involving the links 
Zi... Ik- The links form a pattern which can be as¬ 
signed to a cluster consisting of these links and ad¬ 
jacent sites. If the fluctuation pattern corresponds 
to a disconneted cluster, the overall contribution of 
all summands associated with this pattern annihi¬ 
late each other and such fluctuation patterns can be 
discarded, see also Fig. [l] 

If the links (and the adjacent sites) form a con¬ 
nected cluster ( 7 , the process is called linked and 
the Hamiltonian can be rewritten as sum over these 
linked clusters 


n eS ({*j}) = no+ Y A ? 1 --- A k 

E 7 - kj=k 


v/ci 


. k]\ 


y^C(m) 


/ y —> Y T ™ k ,li ■ ■ ■Tm k ,l k - (20) 

12211 =kYi m i=Q C k h---lk,h\Jh---lk=Ck 


Note that there are also fluctuation patterns which 
correspond to a lower particle channel since not all 
excitations are involved in the fluctuation as illus¬ 
trated in Fig. [3] These fluctuations cancel via the 
subtraction scheme required for the normal order¬ 
ing. 

A linked cluster expansion therefore naturally 
arises from the pCUT approach. The simplest ap¬ 
proach is now to design a sufficiently large cluster 
which incorporates all the relevant linked fluctua¬ 
tions of a desired order and to carry out the calcula¬ 
tion on this single cluster. Yet, the calculation based 
on this approach still incorporates a lot of unlinked 
processes resulting in an unnecessary computational 
overhead. Furthermore, these processes are kept in 
the memory during the calculation which usually de¬ 
fines the highest possible order. 

It is therefore typically favourable to use a full 
cluster decomposition. To this end the contribu¬ 


r 


tion of (small) clusters is calculated by applying H e ff 
to all clusters determining the corresponding matrix 
elements which requires a minimal memory usage. 
Contributions of smaller subclusters have then to be 
subtracted in order to gain the reduced contributions 
of the cluster and to avoid double counting of fluc¬ 
tuations. Hence, after the subtraction, every link lj 
in Eq. (20) is involved at least once. 


The matrix elements in the thermodynamic limit 
are then determined by summing up all the reduced 
contributions of these subclusters. The efficiency of 
this approach is based on the identification of clus¬ 
ters by means of the topology which implies that the 
calculation can be restricted to a much smaller set 
of toplogically distinct clusters called graphs. 
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n x n adjacency matrix mij with 


©= 



Q 




FIG. 2. The fluctuation pattern indicated by the high¬ 
lighted links corresponds to fluctuations associated with 
the vacuum. These fluctuations are irrelevant for the 
local one-particle hopping element in the thermody¬ 
namic limit (left panel). 
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(1,1,0,0,1,0,0,0,1,1)= 803 


FIG. 3. Three different labelings of the same graph are 
depicted with the corresponding adjacent matrices and 
the resulting keys. The first labeling corresponds to the 
canonical labeling and the resulting key is the graph key. 


1. Generating relevant graphs 


In the following, we briefly outline the underlying 
concepts of the toplogical identification by means of 
a canonical labeling. The identification process is 
oriented on the details given in Ref. |T3, 


I 1, if sites i and j are connected by a bond 
mu = < 

|0, else. 

( 21 ) 

Two clusters are called topologically equivalent, if 
a simple renumeration of the sites yields the same 
adjacency matrix in both cases. Obviously, it is suf¬ 
ficient to restrict the calculation to topologically dis¬ 
tinct clusters (graphs) because the results can easily 
be mapped via a relabeling. 

In order to identify clusters as topologically equiv¬ 
alent, it is reasonable to introduce a canonical la¬ 
beling which is distinct from other labelings. To 
this end, one identifies the off-diagonal elements rriij 
(i < j) as the binary bits of an integer number with 
the order (mi^, mi^,... mi ?n ... m n _i ?n ) serving as 
a key. The labeling which maximizes this key cor¬ 
responds to the canoncial labeling and defines the 
graph, see also Fig. [3] for examples. If two clusters 
are topologically equivalent, they are assigned to the 
same graph. In practice, it is necessary to find the 
labeling which maximizes the key efficiently. We do 
not go into technical details for the implementation. 
For further reading we refer to Ref. [13, 

Using the graph key, it is possible to schematically 
produce the required connected graphs by succes¬ 
sively adding links to a given set of graphs. Further¬ 
more, the number of graphs relevant for the calcu¬ 
lation can be additionally reduced. Obviously, it is 
sufficient to include only graphs that can be embed¬ 
ded into the lattice under consideration, i.e. if the 
lattice contains clusters corresponding to this graph. 

By construction, the reduced contribution of a 
graph of n links starts at most at order n. Hence, for 
a given order n, all connected graphs of up to n links 
ahave to be considered. Depending on the problem, 
further selection rules can help to reduce the number 
of relevant graphs because the presence of conserved 
quantities might result in further constraints for a 
given order. 


2. Calculation on graphs 


First, we consider a system where all (super) sites 
are coupled by identical undirected bonds. In this 
case, a cluster is fully determined by an (arbitrary) 
numeration of the sites and the bonds linking the 
sites. A cluster of n sites can be represented by a 


The next step is to perform the actual calcula¬ 
tion on graphs, i.e. H e ^ is evaluated to determine 
the relevant low-energy matrix elements. For the 
implementation, one typically chooses the eigenba- 
sis of Q, i.e. every basis state can be identified with 
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a defined distribution of bare excitations. During 
the application of the operator sequences in H e q the 
intermediate states are represented as linear combi¬ 
nations of these basis states 

i*) = E a iW • ( 22 ) 

3 


The basis states | j) are represented by integer num¬ 
bers where the information of the basis states is en¬ 
coded bitwise in the bit representation of an integer 
as known for instance from exact diagonalization. In 
this basis, the action of the operators in Eq. (18) af¬ 
fects only a small number of bits allowing for a fast 
access and modifications via bitwise operations. 


Compared to the dimension of the full Hilbert 
space, the number of terms in an intermediate state 
is small. Therefore, it is useful to use associative 
arrays consisting of a key (the state represented by 
an integer) and an associated value (the coefficient 
a.j represented by a rational number). This is also 
illustrated in Fig. [4^i. 

The desired matrix elements are then given as se¬ 
ries expansions with rational coefficients. In order to 
gain cluster additive quantities for the embedding, 
the normal ordering is applied followed by the sub¬ 
traction of subclusters contributions. 


IV. WHITE GRAPH EXPANSION 


Up to now, we have described the essential steps 
to set up a linked cluster expansion with the pCUT 
method. This standard scheme becomes inefficient 
when the physical system possesses many linktypes 
corresponding to perturbation parameters A j , which 
is for example relevant when fitting experimental 
data with microscopic calculations. These linktypes 
function as another topological attribute which can 
be incorporated in the usual scheme by generalizing 
rriij G {0,1... n} for n linktypes. Clearly, this leads 
to an exponential growth of relevant graphs needed 
for the calculation. However, as we demonstrate in 
the following, one can circumvent this inconvenience 
by a so-called white-graph expansion. 

The underlying principle of the white-graph ex¬ 
pansion is to ignore the different linktypes (colors) 
for the topological classification of the clusters. The 
white-graph number is simply defined using Eq. © 
and the calculation is restricted to topologically 
distinct clusters in this sense. However, the contri¬ 


butions of clusters in Eq. (20) cannot be gained by 


a simple mapping of the sites because the contri¬ 
butions do differ depending on the color pattern of 
the links. To this end, additional information are 
tracked during the pCUT calculation on the single 
graphs allowing to restore the contribution of each 
cluster exactly from the calculation on white graphs. 


3. Graph embedding 


Finally, the graphs are embedded into the infinite 
lattice, i.e. all realizations of a given graph on the lat¬ 


tice are produced recovering all the clusters in Eq. 20 


via the appropriate relabeling. The contributions of 
the cluster additive quantities are added up accord¬ 
ingly. Note that the contributions of all linked clus¬ 
ters in Eq. (20) are incorporated while the actual 
calculation is performed on a massively smaller sub¬ 
set of topologically distinct clusters. The embedding 
procedure is straightforward and technical details of 
the implementation can be found in Ref. H3l 


Typically, the bottleneck of the calculation is de¬ 
fined by the derivation of the effective model on the 
graphs via pCUTs. Hence, the efficiency of this ap¬ 
proach is based on the identifcation of clusters and 
the associated reduction of the number of clusters 
involved in the calculations. 


A. The calculation on white graphs 

The additional bookkeeping can be interpreted as 
an enhanced or generalized monom; additional val¬ 
ues n are assigned to the links to keep track of the 
relevant information. The latter can be comprised 
in a single object M(n ) which can be viewed as a 
generalized monom. The choice of M(n ) is problem 
dependent and must be adjusted for each applica¬ 
tion. If the links differ only with respect to the cou¬ 
pling strength, it is sufficient to track how often each 
link was active in the operator sequence applied to 
the current state. For more complex problems, more 
information must be tracked and this approach may 
become unfavorable in scenarios where a lot of link- 
types are involved and the non-zero matrix elements 
differ greatly between the linktypes. However, in 
principle it is always possible to track the relevant 
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(a) 



+ 



+ ... 


FIG. 4. (a) The stan dar d representation of intermediate 
states following Eq. ( [22] ) so that ctj corresponds to a ra¬ 
tional number, (b) Standard calculation on white graps 
expansion using the generalized monoms as written in 
Eq. ( [24| . (c) Improved calculation on white graphs fol¬ 
lowing Eq. (25) which uses the product structure of the 
generalized monoms. 


information. 

Evaluated for a given link pattern, the generalized 
monoms must yield an actual polynom in the mul¬ 
tiple perturbation parameters of the different cou¬ 
plings such that the result corresponds exactly to the 
calculation on a cluster having this very link pattern. 
This allows for the subsequent evaluation necessary 
for the embedding procedure. 

As indicated by the name, the generalized monoms 
obey 


M(n 1 ) • M(tT 2 ) = M(n 1 + n 2 ) (23) 

and the intermediate states during the calculation 
are then expressed as linear combinations of these 
elements: 


states. The information of a generalized monom can 
analogously be encoded bitwise. In a naive imple¬ 
mentation, the state combined with the generalized 
monom defines the key and the associated amplitude 
aij corresponds to the value as visualized in Fig. EF 
In an improved implementation, one can make use 
of the product structure of the generalized monoms 
as well as of the intermediate states and factorise the 
representation: 


W = XT’) ( 


(25) 


The effect of the operators on the state | j) can be 
calculated independendly of the effect on the appen¬ 
dant sum of generalized monoms. Consequently, nu¬ 
merous redundant search, comparison, and shift op¬ 
erations are avoided. 

The implementation of this approach requires two 
nested containers. The key of the first container is 
defined by the state | j) while the value is given by a 
second container. The second container is defined by 
a key comprising the information of the generalized 
monom M(n$) and the value corresponds to the am¬ 
plitude aij of this generalized monom and the state. 
This is illustrated in Fig. EF 

For an efficient performance, we recommend a sim¬ 
ple modification of this procedure. As shown in sub¬ 


section mm, the reduced contribution of a graph 
involves each link at least once. This property can 
help to reduce the computational effort immensely. 
The objective is to calculate directly the reduced 
contribution of a cluster without relying on a sec¬ 
ond subtraction step. 

We consider the calculation in order k on a graph 
consisting of n links. Let denote how often a 
link j appeared in the operator sequences comprised 
in the monom M(n$). These information are typ¬ 
ically tracked anyway and cause no computational 
overhead. We define 




(24) 




Note that the action of an operator in Eq. (18) af¬ 


fects the generalized monoms M(f^) and the state 
| j) independently; similarly to the action on product 
wave functions. 

As in a standard pCUT set up, it is reasonable to 
use associative arrays to represent the intermediate 


(Jii = ^ fiij with 
j 


j^i,j — /Al 
k*i,j — 0 , 


if Hij > 1 
else. 


(26) 


One can discard all intermediate generalized 
monoms M(n^) where fii > k — n. 

Evidently, this is specifically relevant if the number 
of links n is close or equal to the order fc, i.e. for 







9 


graphs with a large number of links. By construc¬ 
tion, these graphs represent the majority of graphs 
making this modication extremely valueable. 

Finally, the matrix elements of the Hamiltonian 
are given as linear combinations of the generalized 
monoms. The generalized monoms of these matrix 
elements build the centerpiece of this approach car¬ 
rying the relevant information for the embedding of 
the white graphs described in the following. 


B. Embedding white graphs 


For the embedding procedure of a standard linked 
cluster expansion, all realizations of a given graph on 
the lattice are produced recovering all the clusters in 
Eq. (20) via the according relabeling and the contri¬ 


butions of the clusters are added up. Analogously, 
for the white-graph embedding, all realizations of a 
given graph on the lattice are generated. In order to 
recover the contributions of all clusters in Eq. [20| the 
information of the link pattern must be incorporated 
in a second evaluation step. The matrix elements on 
the white graphs are given as sums of generalized 
monoms. Evaluated for a given link pattern, the 
generalized monoms yield an polynom in the mul¬ 
tiple perturbation parameters recovering the result 
of the calculation on a cluster corresponding to the 
given link pattern. The resulting contributions are 
added up for all embeddings recovering the correct 
result in the thermodynamic limit. 


V. COUPLED TWO-LEG XXZ 
HEISENBERG LADDERS 


In order to illustrate the functioning and useful¬ 
ness of a white-graph expansion, we discuss in this 
section the calculation of one-magnon dispersions for 
the ordered state of a two-dimensional quantum spin 
model of coupled two-leg XXZ ladders involving four 
different pertubation parameters. The microscopic 
model in terms of spin 1/2 operators reads 


H = 


X 

7, (i,j) 


SfS! + ±(S?Sf + S?S?) 


• (27) 


and is illustrated in Fig. [5] The different couplings J 7 
correspond to rung (legjexchange J rung ( J\ eg ) of the 
two-leg ladders, the inter-ladder exchange Ji nt , and 



FIG. 5. Illustration of the quantum spin model of cou¬ 
pled two-leg XXZ ladders. Filled circles represent quan¬ 
tum spin 1/2 which are linked by various solid lines illus¬ 
trating the different exchange couplings Jj_, Jy, and Ji n t- 
Black box shows the two-site unit cell of the system cor¬ 
responding to the rungs of the ladder and 5x ( Sy ) refers 
to the distance between two neighboring rungs belonging 
to different (same) ladders. 


a spin-anisotropy A. This Hamiltonian was recently 
shown to be relevant for the experimental quantum 
magnet CgHig^CuB^ displaying long-range Neel 
order and gapped magnon excitations m- 

To calculate the spin-wave dispersion, we perform 
a sublattice rotation to obtain a ferromagnetic refer¬ 
ence state for the Ising case A = 0, which transforms 
Eq. ( |27| ) to 


7 = X J 7 

7, (hi) 


_ QZ Q Z 

^ j 


(S+S+ + Sr S r) 


(28) 

We then introduce hardcore-boson operators a\ and 
a v (bl and b u ), which create and annihilate a magnon 
at site A ( B ) of rung v above the ferromagnetic ref¬ 
erence state, obtaining the hardcore-boson Hamilto¬ 
nian 



+ 53 [T„,o + A (T„,_ 2 + T„ m )\ , (29) 

V 


where J = J\ eg + ( Jrung + Ant)/2, N is the number of 
unit cells, n^ = a\a v , and nfi = b\b v . The sums 
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are taken over all rungs. The operators T^ n , with 
T v - 2 = 2 , are given by 

T _ — _ r 

- L v , ,0 ^rung'^v ^lnt'^v L v-\-8x 

-*leg (nl^nl+Sy + ( 30 ) 

and 

^k,+2 = ^rungh^S^, ^int^t^I+(5cc 

-*ieg » ( 31 ) 


where x 1 = J 7 /j, and &e (5i/) is the distance be¬ 
tween two neighboring rungs belonging to different 
ladders (the same ladder). 

At this point let us replace x 1 —>• rx 1 with 
r G [0,1] so that r = 1 corresponds to our physical 
system Eq. (29) and r = 0 functions as a well defined 


starting point for perturbation theory. Indeed, the 
resulting Hamiltonian can be rephrased as 


J = n 0 + r £ : 

J ne{ —2,0,2} 


(32) 


meeting all criteria relevant to apply the pCUT 
method as described in detail above: i) The unper¬ 
turbed Hamiltonian at r = 0 is equidistant, ii) the 
unperturbed spectrum is bounded from below, and 
iii) the perturbation is decomposed in T n operators. 
One can therefore map Eq. ( |32| to an effective model, 
% e ff, which conserves the number of magnons. 

The one-magnon sector 'H^, which is of our in¬ 
terest here, corresponds in real space to an effective 
hopping Hamiltonian for the magnons and is there¬ 
fore fully determined by the one-magnon hopping 
amplitudes. Here we have calculated all these hop¬ 
ping amplitudes by series expansion up to 13th order 
in all parameters rr 7 . At this order one has to treat 
2709 white graphs in total. Clearly, the number of 
graphs with color is many orders of magnitude larger 
and a related calculation is far beyond any realistic 
set up. 

Since the different linktypes differ only with re¬ 
spect to their coupling strengths, the additional in¬ 
formation n tracked during the calculation on white 
graphs must incorparate how many times each link 
is acted on by the perturbation parameters J 7 . This 
is in fact identical to the approach where all links 
of a graph have different coupling constants and the 
result is given as a multi-variable polynom. Finally, 











J\\JL 


h 

h 

h 


FIG. 6. Illustration of two different fluctuation patterns 
contributing to the local one-particle hopping element 
tij in thermodynamic limit (left panel) for the quantum 
spin model of coupled two-leg XXZ ladders. The solid 
lines correspond to the different exchange couplings J± 
(red), J|| (b lue ), and Jmt (green). The occurrence of a 
link in Eq. ( |20| ) is visualized by an arc defining the fluc¬ 
tuation pattern. All fluctuations with this pattern are 
pooled together in this representation. Both fluctuations 
are associated with topologically distinct graphs Gi and 
C /2 (right panel) solely because their colorings are dif¬ 
ferent. Decomposing the results of the calculation with 
respect to the fluctuation pattern allows to gain both 
contributions by embedding a single white graph C/ whlte 
and ’evaluating’ the resulting coloring. The evaluation 
yields polynoms in the different perturbation parameters 


during the embedding procedure, the coupling con¬ 
stants of the graph are matched with the ones of the 
specific graph realizations on the lattice. This is vi¬ 
sualized for the local one-particle hopping element 
ti,i in Fig. [6j 

The one-magnon sector can then be simplified 
by Fourier transform to 

^eff = 5-/ (^ a ^ a k a k + 5 (33) 

k 

where uj a (k) and uj^{k) denote the two one-magnon 
branches stemming from the two-site unit cell. 

Let us stress that the results for the one-magnon 
dispersions are given analytically in all physical pa¬ 
rameters x±, £||, Xi n t, and A. Scanning through a 
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large variety of parameter sets, which is often needed 
when fitting experimental data, is then straigh- 
foward as long as convergence in r up to r = 1 is 
given: the values for the different couplings just have 
to be inserted into the analytical expression. 

Physically, the series in r do converge well for A 
small. First, both limits r = 0 and A = 0 have 
exactly the same ground state as well as the same 
one-magnon excitations. Only the multi-magnon 
energies differ. Second, quantum fluctuations are 
strongly suppressed for small A and/or r. We find 
that bare series in r are quantitatively converged up 
to r = 1 for A < 0.8. 

For larger values of A, one has to rely on extrapola¬ 
tion schemes like dlogPade extrapolation. This is ex- 
pecially true for the most demanding case of the one- 

magnon low-energy gap A = min^ (cu a (k), cu^(k)^ at 

A = 1. In dlogPade extrapolation one constructs var¬ 
ious extrapolants [L, M] where L ( M ) denotes the 
order of the numerator (denominator). Explicitly, 
the dlogPade extrapolation is based on the Pade ex¬ 
trapolation of the logarithmic derivative of the one- 
particle gap A 


d 

— In A 
dr 


[L,M] 


Pl_ 

Qm 


(34) 


where Pl and Qm are polynomials of order L and M. 
Due to the derivative of the numerator in Eq. [34] one 


requires L + M = m — 1 where rri denotes the max¬ 
imum perturbative order which has been calculated 
(here m = 13). The [L, M] dlogPade extrapolant is 
then given by 


iL ' Mi:=exp (i w) dT ) - <35> 

In the case of a physical pole at To = 1 one is 
able to determine the dominant power-law behaviour 
\t — ro\ ZL/ close to To- The exponent zv is then given 
by the residuum of Pl / Qm at t = tq 


zv = 


Pl(t) 


T = Tq 


(36) 


The quality of the dlogPade extrapolation can be 
further improved when the location of a quantum 
critical point is known exactly. This is true in our 
case since the only physical point of the expansion 
is tq = 1. For this case one can bias the dlogPade 



FIG. 7. One-magnon gap A/J as a function of r for 
the isotropic square lattice J\ eg = J rung = Jint = A = 1 
for different bare series and different (biased) dlogPade 
extrapolants. Inset Zoom close to r = 1. 


extrapolation by setting 


B [L. M] := exp 



W) i 

Qm(t') (r - To) 



where 


(37) 


(r-T 0 ) s !nA 


[L,M] 


Pl_ 

Qm 


(38) 


All bias dlogPade extrapolants B[L,M] have a pole 
at r = ro by construction. The value of the denomi¬ 
nator of B [L, M\ at r = ro corresponds then to the 
critical exponent zv. 

In the following we exemplify a typical behaviour 
of A for the isotropic Heisenberg model on the square 
lattice with J\ eg = J rung = J- m t — A = 1. This 
isotropic point is known to have gapless Goldstone 
bosons due to the breaking of SU(2) symmetry in the 
long-range ordered Neel ground state. One therefore 
expects that the one-magnon gap A is gapped for all 
values r < 1 and vanishes at r = 1. As outlined in 
Ref. 0 the quantum critical behavior is mean-field 
like, i.e. one has critical exponents v — 1/2 and 2 = 
1. Our results for the one-magnon gap are displayed 
in Fig. [7] The bare series is reliable up to r < 0.85 
predicting an unphysical gap for r = 1. 

This situation is strongly improved using dlogPade 
extrapolation. Essentially all extrapolants show a 
pole close to r « 1. Taking all [L, M\ of highest or¬ 
der L + M = 12 with L > 3 and L > 3, the average 
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pole is found at To = 1.004 with an average critical 
exponent zv « 0.52 very close to the expected value 
1/2. The quality of the extrapolation gets even bet¬ 
ter when biasing the extrapolants such that they ex¬ 
hibit a pole at exactly r — 1. In this case one finds 
an average critical exponent zv ~ 0.499. Overall, 
we therefore find quantitative agreement even in the 
most demanding quantum critical regime using our 
high-order series expansion. Let us stress that simi¬ 
larly well converged results can be directly deduced 
in all parts of the (large) parameter space. 

VI. CONCLUSIONS 

In this work we have introduced the so-called 
white graph expansion for the method of perturba¬ 
tive continuous unitary transformations when imple¬ 
mented as a linked cluster expansion. The essential 
idea behind this expansion is an optimal bookkeep¬ 
ing during the calculation on graphs which is possible 
due to the presence of the model-independent effec¬ 
tive pCUT Hamiltonian in second quantization. Our 
approach is especially useful when the microscopic 
model under consideration has several expansion pa¬ 
rameters. This case usually represents a major chal¬ 


lenge for any kind of linked cluster expansion due to 
the proliferating number of linked graphs. Indeed, 
each expansion parameter corresponds to a distinct 
color in the graph expansion. The white graph ex¬ 
pansion overcomes this challenge to a great extent, 
since the actual calculation is performed on white 
graphs and the coloring is done after the calculation 
as a final step. 

We are strongly convinced that this white graph 
expansion is useful in many situations, e.g. micro¬ 
scopic models in three dimension or systems with 
long-rang interactions being notoriously complicated 
for linked-cluster expansions. Finally, it would be 
interesting whether other types of linked-cluster ex¬ 
pansions like for expample Takahashi’s degenerate 
perturbation theory |38| or matrix perturbation the¬ 
ory [13] can be also reformulated as a white graph 
expansion. 
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